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Efficient Calculation of Current Flow in 
Electroplating Cells 
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(Manuscript received August 31 , 1979) 

An efficient computational method has been developed for calcu- 
lating the current density distribution within electroplating cells. The 
computer program can analyze two-dimensional cells with arbitrary 
polygonal shapes, as well as three-dimensional axisymmetric cells 
with arbitrary polygonal cross section. Laplace's equation in the bulk 
of the cell is coupled to the convectiue diffusion equation in the 
boundary layers near the electrodes. Finite-difference and finite- 
element methods, and infinite sums, are avoided by a boundary 
integral formulation for the solution of Laplace's equation. The 
method is illustrated with sample calculations for copper plating of 
multilayer boards. The predicted copper plating thickness distribu- 
tion is essentially the same as the primary current distribution. 

I. INTRODUCTION 

Until recently, analytic and numerical techniques have not been 
available for analyzing electroplating cells. Even after the work of 
Newman and coworkers, 1 only restricted geometries were practical. 
Design of new electroplating systems has therefore necessarily been 
done by educated guess, rule of thumb, extrapolation from previous 
systems, and expensive experimentation. Although the physical and 
chemical processes taking place in a real-life electroplating cell are 
generously complex and nonlinear, realistic mathematical models can 
be constructed, and it is within the state of the computing art to solve 
these models efficiently. 

This paper describes a computer program to solve the mathematical 
models describing a wide variety of electroplating cells. An important 
consideration during this development has been that the resulting 
computer program be able to analyze many different cell geometries 
and flow patterns with at most trivial program changes. A second 
consideration has been that the program be efficient and inexpensive 
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to run, since many geometry and parameter studies are likely to be 
made in the design or analysis of an electroplating system. 

The method is illustrated with sample calculations for copper plating 
of multilayer boards. 

II. PROBLEM FORMULATION 

A schematic cross section of a typical copper-plating tank is shown 
in Fig. 1. Copper anodes both provide current flow and replenish the 
plated metal in the solution. Cathodes, the multilayer boards, receive 
the plated metal and provide the return path for the current. The 
electroplating solution contains water, a salt of the plating metal, in 
this case CuS0 4 , an acid to increase the conductivity of the solution, 
in this case H 2 S0 4 , and proprietary secret ingredients which improve 
the microscopic structure of the plated metal. The solution is well- 
stirred, and a flow of fresh solution past the cathodes is provided. In 
the typical tank, there are three different stirring mechanisms. First, 
the cathodes may be rocked back and forth along their width (perpen- 
dicular to the plane of Fig. 1) to reduce nonuniformity of plating from 
the bar anodes. Second, solution is pumped out of the tank near the 
bottom, filtered to remove foreign matter and granular copper, and 
then returned to the tank. Finally, each cathode has a sparger, or air 
bubbling pipe, under it, directing a stream of air bubbles up each side 
of the cathode and producing the convection pattern sketched in 
Fig. 1. 

Other complications are, of course, possible and may be analyzed by 
the same program. Anodes and cathodes may be of more complex 
shapes; other metals or insulators may also be present, to alter the 
fluid flow or the distribution of plated metal on the cathode. 
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Fig. 1— Schematic cross section of a multilayer board plating tank, showing anodes 
(A), cathodes (C), and spargers (S). 
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Near the cathodes, bubbles from the spargers rise, causing a convec- 
tion cell between each cathode-anode pair. For analyzing many of the 
possible effects, such as geometry, it is sufficient to assume a periodic 
tank, and to analyze only one convection cell, as in Fig. 2. (The 
program can analyze two-dimensional cells with arbitrary polygonal 
shapes, as well as three-dimensional axisymmetric cells with arbitrary 
polygonal cross section.) Very thin boundary layers, not shown in the 
figure, are adjacent to the electrodes. In the boundary layers, the 
copper concentration obeys a diffusion-convection equation. For ve- 
locity flows of the boundary layer type, 



V X = 2YE(X) 



y Y __y 2 dE(X) 



dX 



(1) 
(2) 



the diffusion-convection equation in the boundary layer can be inte- 
grated approximately, 2 giving the concentration in terms of the current 
flow at the electrode. 

The velocity fields in the boundary layers are needed. The velocity 
flow caused by bubbles rising near a vertical flat plate is complicated, 
and the solution for the flow is not known. However, the flow is 
qualitatively similar to the case of uniform flow past a flat plate, whose 
solution is available (Ref. 3, p. 125). If the flow far from the plate is V. 
in the Y direction, the velocity field is as above, with 

E(X) * 0.166[V 3 „/(vA")] 1/2 . (3) 



j = 



j = 



i = o 



j = 



Fig. 2— One electroplating "cell." The parts of the boundary with no current flow are 
marked by j = 0. 



CURRENT FLOW IN ELECTROPLATING CELLS 171 



With s measuring distance along the cathode, the methods of Appendix 
II of Ref. 3 give 



-if 



j(t) dt 

(s a7r =Wy 



C ' S ' ' i I z-,3/4 _ *3/4\2/3 W 



with the dimensionless parameter k c given by 

'32D 2 0.166 V^T 73 T(%)C b F 2 XV 2 



k c = 



kR b T 



(5) 



g- 



The velocity field at the anodes is more difficult, since each anode 
is enclosed in a porous bag. The flow must be some kind of free 
convection, caused by the variation of electrolyte density with copper 
concentration. Since the only item of interest is the copper deposition 
distribution on the cathode, it seems unnecessary to include all the 
complications at the anode. Rather than assume a velocity field, a 
simple model for the concentration has been adopted. If s measures 
distance along the anode, from the bottom, 

c(«)-l+7^j. (6) 

k a S 

This approximation is equivalent to assuming that the anode boundary 
layer has width proportional to s* and that the concentration varies 
linearly across the boundary layer. The effect on the cathode current 
distribution of various choices of k„ and 8 can be investigated later. 
Reasonable values of S probably are in the range of from Vt (free 
convection) to Vfc (uniform flow past a flat plate). 

The derivation of the equations describing copper plating is by now 
standard 1,2 and will not be repeated. The dimensionless variables we 
will need are the potential, <J>, the concentration of copper, c, and the 
current density flowing into the cathode, j. In our normalization, j is 
equal to the directional derivative, 

>-£■ 

where n is the outward-pointing unit normal vector at each point of 
the cathode. The normalizations used are Xo for lengths, R M T/F for 
potentials, kR^T/XoF for currents, and C/,, the bulk concentration, for 
the copper concentration. 

Throughout the region shown in Fig. 2, outside the boundary layer, 
Laplace's equation holds: 

V 2 = 0. (8) 

On the anode, the potential is a given constant, <t> = <j> a . The cathode 
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is grounded. On the remaining boundaries, no current flows, d$/dn = 
0. All the complication is in the nonlinear boundary conditions on the 
electrodes, coupling <£ and j, in effect defining a nonlinear resistance of 
each boundary layer. The potential across the cathode boundary layer 
has two constituents, 

4> = tj s + TJc. (9) 

For the copper electroplating problem of interest, the concentration 
overpotential is 

j] c = '/ 2 ln(c). (10) 

The surface overpotential is given by the relation 

j = k e c y (e a " n * - e-"**) , (11) 

with the dimensionless parameter k e given by 

_ FXolo 
ke kR h T ■ (12) 

The boundary layer at the anode is treated similarly, although in 
practice the details are not important. 

The problem is thus reduced to solving Laplace's equation in a 
rectangular region. There are linear boundary conditions on part of 
the boundary and complicated nonlinear boundary conditions on the 
remainder. An efficient technique for solving this problem is outlined 
in the next sections. 

III. SOLUTION OF LAPLACE'S EQUATION WITH NONLINEAR BOUNDARY 
CONDITIONS 

We wish to solve Laplace's equation on a two-dimensional region, 
with boundary B. Let s denote arc length on B. On the electrodes, a 
nonlinear boundary condition couples <f>(s) and j(s) = d<f>/dn. On the 
rest of the boundary, j{s) = 0. 

We first discuss solving Laplace's equation with nonlinear boundary 
conditions, assuming solutions with linear boundary conditions are 
available. In the next section, Laplace's equation with linear boundary 
conditions is discussed. 

We expand the potential in terms of known solutions, 

/, 

<t> = S di4> h (13) 

/=i 

and determine the L coefficients d/ so that the nonlinear boundary 
conditions are satisfied. Letting j) = d(f)//dn, <£>/ is defined to be a 
solution of Laplace's equation with <j>/ = b/(s) on the electrodes, and 
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ji(s) = on the rest of the boundary. The functions bi(s) are in principle 
arbitrary; fl-splines are used as discussed in the next section. 

Also choose M points s m on the electrodes. M > L is required; we 
use M = 2L. Suppose approximate values for the rf/s are known. Then 
j(s) = ££-i dji(s) is the approximate current. From the integral formula, 
calculate the concentration c(s m ) at M points; then calculate the 
concentration overpotential -r\ c (s m ). Then obtain the surface overpoten- 
tial Tj s (s m ) by a Newton iteration. Finally, at M points do a least-squares 
fit to 

rjc(s) + i,.(8) = J dibds). (14) 

The differences, di — di, are zero when the correct boundary conditions 
are obeyed. 

The preceding paragraph defines a set of L nonlinear equations for 
the d"s. An initial guess is made, and the equations solved by a quasi- 
Newton iteration with rank-one update. 4 ' 5 Usually only a few iterations, 
typically 5 to 10, are needed. Convergence is somewhat slower when 
the current is large, since the equations then are effectively more 
nonlinear. 

In practice, the nonlinear equations are considerably less sensitive 
to the values of the di's on the anode than to those on the cathode, 
especially for high currents. This leads to difficulties in solving, since 
the Jacobian matrix is therefore ill-conditioned. This has been over- 
come by a two-step iteration process. First, <f> on the anode is held 
fixed, and Laplace's equation is solved with nonlinear boundary con- 
ditions on the cathode. Second, the resulting <f> on the cathode is held 
fixed, and Laplace's equation is solved with nonlinear boundary con- 
ditions on the anode. These two steps are repeated. Only two or three 
of these two-step iterations are needed. 

IV. LAPLACE'S EQUATION WITH LINEAR BOUNDARY CONDITIONS 

To be able to deal with a wide variety of geometries, a boundary 
integral method is used. A more detailed discussion appears in Ref. 6. 

We wish to solve V 2 <J> = in a two-dimensional region with boundary 
B. On part of B, B\, the boundary condition is <f> = u(s); on the 
remainder, B 2 ,j= d<f>/dn = 0. 

We start with Green's third boundary identity 7 



2™j>(x, y) = 



8C 

<i>(s) — -j(s)G 
an 



ds. (15) 



Here (x, v) is some point inside B, and G is the "fundamental solution" 
to Laplace's equation, 

G = -Vi Hi** ~ x) 2 + (y s - y) 2 ]. (16) 
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Also, (x Hi y„) are the coordinates of the point at arc length s on B, and 
n(s) is the outward normal at s. In three dimensions with axial 
symmetry, the formulation is identical except for a different funda- 
mental solution. 

If <f>(s) andy'(s) are replaced by other functions <f)*(s) andy'*(s), the 
function # defined by 



2ir<t>(x, y) = 



dG 

'(s) — -j*(s)G 
an 



ds (17) 



exactly obeys Laplace's equation inside B, but does not obey the 
correct boundary conditions. The linear mixed boundary problem may 
be solved approximately by choosing functions <f>*(s) andy'*(s) which 
make </> approximately obey the boundary conditions. 

To do this, we parameterize <f>* and 7* and determine reasonable 
values for the parameters. On B } we take 

n, 
**(*)- u(s) and /*(«) - £ Okbkis). (18) 

On B 2 we take 

N t +N 2 

<(>*(s)= X akbk(s) and j*(s) = 0. (19) 

The a's are coefficients to be determined so that the boundary condi- 
tions are satisfied as closely as possible. The b k may be any useful set 
of functions. We use B-splines. 8,9 As a basis function for./*, a function 
behaving like the inverse square root of distance from the lower end of 
each electrode may be used, in addition. Such a term is present in the 
primary current, though not in the solution to the nonlinear problem. 
Inclusion of this basis function affects the result only very near the 
end of the electrode; its inclusion assures that far fewer basis functions 
are necessary for good accuracy. This basis function is included in the 
sample calculations; its use in the program is optional. 

If (x, y) approaches, from inside B, a point at arc length t at a 
smooth part of B, it may be shown 7 that the 2tt in (17) changes to a m 
and the line integral becomes a Cauchy principal value integral. Then 
we have 



t on B\:iru{t) 

N > +N * J v a / . dG ^ 

*onB 2 :77 X aMt) "J ^ akbk{s) ^ ds 



ton Bi'.mi(t) 

I 

N,+N 2 ~| 

ds, (20) 
where the integrals are to be interpreted as Cauchy principal value 



+ 



dG N,+N 2 

u(s) — -G 2 a k b k (s) 
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integrals at s = t. For any particular point tj, this equation reduces to 
the form 

N,+N 2 

r,= £ Ajkdk, (21) 

where Ajk involves integrals of b k (s), and r, an integral of u(s). Some of 
the integrals are singular, and require some care. 6 Equation (20) cannot 
be made to hold for all t on B, since there are only Ni + N 2 a's. An 
approximate solution for the a's can be found by choosing Na points tj, 
Na> Ni + N 2 , and doing a least-squares fit of the Afc linear equations 
in N\ + N2 unknowns. 

The resulting a*'s give <f>*(s) and j*{s) directly, with no need to 
perform any derivatives. An error estimate for may be obtained by 
comparing <f>*(0 and the <f>(t) obtained from the integral definition, 
letting t -> B. 

In practice, a second integral equation is also used, to improve the 
conditioning of the A jk matrix. Details are in Ref. 6. 

V. SAMPLE CALCULATIONS 

In this section, sample calculations are presented for a typical tank. 
Its dimensions are cathode length 46 cm, anode length 46 cm, cathode- 
anode separation 16 cm, and total tank depth 84 cm. We use X = 46 
cm. The three major components of the electrolyte are Cu ++ , H + , and 
HSCV The bulk concentrations are approximately 0.27 molar CuS0 4 
and 1.7 molar H 2 S0 4 . The bulk solution has a concentration of copper 
ions of C h = 0.00027 moles/cm 3 . The bulk conductivity, k, is obtained 
using data in Newman 1 for mobilities, giving k = 0.63(ohm — cm) -1 . 
Also D = 0.55 X 10~ 5 cm 2 /s and v = 0.01 X 10" r> cm'Vs. Then the 
dimensionless parameters are a, ft, y, S, k a , k c , and k e . The remaining 
material constants are taken to be: 10 a c = 1.5, a a = 0.5, and y = 0.42. 
Also from Ref. 10, I = 0.001 amp/cm 2 . We obtain k c s 26 Vl /2 and k c 
= 2.5. 

Using k c = 260 and <f>„ = 30, the first set of sample calculations tested 
the dependence of j(x) and c(x) on the anode parameters k a and S. 
With k„ = 50, values of 5 used were -V2, -V*, 0, and +V6; with S = - l A, 
values of k a used were 10, 50, and 250. As hoped and expected, j(x) and 
c(x) at the cathode were quite insensitive to the choices of the anode 
parameters k„ and 5, varying by only a few percent over all the choices 
of anode parameters. For the remainder of the sample calculations, k a 
= 50 and 5= -Va. 

The second set of calculations varied the applied potential, <j> a , 
keeping k c = 260. Figures 3 and 4 show>(x) and c(x), for <f> a = 20, 30, 
and 40, as well as their limiting distributions and the primary current 
distribution. (In Figure 3, the curve for </>„ = 30 has been omitted; it is 
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NORMALIZED DISTANCE ALONG CATHODE 

Fig. 3 — Normalized plating thickness distribution on cathode for anode potentials 
(p„ = 20 and 40. Other parameters are k„ = 50, S = — Vi, and k, = 260. For comparison, the 
primary distribution (P) and the limiting distribution (L) are also shown. 
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Fig. 4 — Copper concentration on cathode for anode potentials <£„ = 20, 30, and 40. 
Other parameters are k„ = 50, S = — '/«, and k, = 260. 
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Fig. 5— Normalized plating thickness distribution on cathode for k, = 170, 260, and 
520. Other parameters are <>„ = 30, k„ = 50, and S = -V*. 
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Fig. 6— Copper concentration on cathode for k, = 170, 260, and 520. Other parameters 
are <£„ = 30, k„ = 50, and S = -'/«. 



almost identical to that for fa = 20.) Each j(x) curve has been 
normalized to its mean value, to emphasize the shapes of the curves. 
The primary current distribution is the current derived from solving 
Laplace's equation with linear boundary conditions, assuming zero 
potential drop across the concentration boundary layers. Recall that 
j(x) is proportional to the rate of deposition of copper on the cathode, 
and therefore is proportional to the plating thickness distribution. As 
the applied potential increases, the concentration of copper in solution 
at the cathode decreases. This decrease does not significantly change 
the shape of the j(x) curve from fa = 20 to fa = 30. At fa = 40, c(x) 
becomes quite small near x = 1; the potential between the bulk and 
the cathode, tj c (x) + tj s (x), becomes large; and they're) curve changes 
shape. With still higher applied potential, j(x) becomes less uniform. 
A limiting current is defined as the j(x) in (4), which makes c(x) 
identically zero, 



lUx)= 4rmr(Vs)x^ x0 - 207 x- 



1/2 • 



(22) 



Since t] c (x) is proportional to In c(x), the limiting current requires 
infinite applied potential. For Fig. 3, the total currents are approxi- 
mately 35, 58, and 79 percent of the total limiting current, for fa = 20, 
30, and 40, respectively. Practical copperplating usually is done at a 
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Fig. 7— Normalized plating thickness distribution on cathode for anode lengths of 41 
cm (S), 46 cm (E), and 51 cm (L). Other parameters are <i>„ = 30, k„ = 50, S = -Va, and 
k c = 260. 
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Fig. 8— Copper concentration on cathode for anode lengths of 41 cm (S), 46 cm (E), 
and 51 cm (L). Other parameters are <t>„ = 30, k„ = 50, S = - l A, and k, = 260. 



total current less than one-quarter of the total limiting current. These 
elaborate calculations demonstrate that, for these typical copperplat- 
ing conditions, the primary current distribution is quite adequate for 
predicting the plating thickness distribution. In the remaining calcu- 
lations, <$> a = 30 is used. 

The third set of calculations varies the cathode parameter k c , which 
is proportional to the square root of the fluid velocity. Figures 5 and 6 
show the results. If k c is too small, the current becomes quite nonuni- 
form. For the remaining calculations, k c = 260. 

The final set of calculations shows the effect of anode length. Anode 
lengths of 41, 46, and 51 cm were used; the cathode length was left at 
46 cm. Figures 7 and 8 show the results. With the shorter anode, the 
current was within ±10 percent over more of the cathode than with a 
longer anode. 



180 THE BELL SYSTEM TECHNICAL JOURNAL, FEBRUARY 1980 



APPENDIX 

Symbols 

at, Unknown coefficient in solution of Laplace's equation. 

Ajk Matrix element. 

bi Function in expansion of boundary values. 

B Boundary of region. 

B\ Part of boundary on which $ is given. 

Bi Part of boundary on which no current flows. 

c Normalized copper concentration, dimensionless. 

Cb Bulk copper concentration, g-mole/cm 3 . 

d/, d'i Coefficients in expansion of overpotential. 

D Diffusion coefficient, cm/s. 

E Velocity field function, s" 1 . 

F Faraday's constant, 96500 coul/g-equivalent. 

G Fundamental solution to Laplace's equation. 

To Exchange current density, A/cm 2 . 

j,j* Normalized current density, dimensionless. 

ji Normalized current density of /th solution, dimension- 
less. 

yiim Limiting current distribution, dimensionless. 

k a ,k c , k c Dimensionless constants. 

L Number of coefficients in expansion of overpotential. 

M Number of fitting points for overpotential expansion. 

N\ Number of coefficients in expansion of j on B\ . 

N2 Number of coefficients in expansion of <f> on Bi . 

N.i Number of fitting points in solution of Laplace's equa- 
tion. 

n Directional derivative, in d<p/dn and dG/dn. 

n Unit vector in outward normal direction. 

Tj Matrix element. 

Rf, Universal gas constant, 8.21 joules/g-mole deg. 

s, t Arc length along boundary, dimensionless. 

s m Fitting points on boundary. 

T Absolute temperature, degrees Kelvin. 

u Boundary condition for <>. 

Vx, Vy Components of fluid velocity, cm/s. 

V x Fluid velocity far from plate, cm/s. 

x, y Coordinates, dimensionless. 

x s , y s Coordinates of point at distance s along electrode. 

X, Y Coordinates, cm. 

Xq Normalizing length, cm. 

ota, etc, y Parameters in reaction rate expression, dimensionless. 

8 Exponent in concentration at anode, dimensionless. 
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r Gamma function. 

7j c Concentration overpotential, dimensionless. 

tj s Surface overpotential, dimensionless. 

k Bulk electrolyte conductivity, (ohm — cm) -1 . 

v Viscosity, cm 2 /s. 

<£, <j>* Normalized potential, dimensionless. 

<£a Anode voltage, dimensionless. 

<j>i Potential for Ith solution of Laplace's equation, dimen- 
sionless. 
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